#include <errno.h>
#include <math.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "bezier.h"

#define ERR(e, r)		\
	do {			\
		errno = (e);	\
		return (r);	\
	} while (0)

static int	addvalues(double **, size_t *, size_t, double *,
		    size_t);
static ssize_t	quadratic(double, double, double, double *, int);

/* Evaluate for X at `t' */
ssize_t
bz1x(double *bz, double *o, double *t, size_t n)
{
	double a;
	double b;
	size_t i;
	ssize_t r = 0;

	if (bz == NULL || t == NULL)
		ERR(EFAULT, -1);
	if (o == NULL) {
		for (i = 0; i < n; i++)
			r += t[i] >= 0 && t[i] <= 1;
		return (r);
	}
	a = bz[2] - bz[0];
	b = bz[0];
	for (i = 0; i < n; i++) {
		if (t[i] < 0 || t[i] > 1)
			continue;
		o[r++] = a * t[i] + b;
	}
	return (r);
}

/* Evaluate for Y at `t' */
ssize_t
bz1y(double *bz, double *o, double *t, size_t n)
{
	double a;
	double b;
	size_t i;
	ssize_t r = 0;

	if (bz == NULL || t == NULL)
		ERR(EFAULT, -1);
	if (o == NULL) {
		for (i = 0; i < n; i++)
			r += t[i] >= 0 && t[i] <= 1;
		return (r);
	}
	a = bz[3] - bz[1];
	b = bz[1];
	for (i = 0; i < n; i++) {
		if (t[i] < 0 || t[i] > 1)
			continue;
		o[r++] = a * t[i] + b;
	}
	return (r);
}

/* Evaluate for X at `t' */
ssize_t
bz2x(double *bz, double *o, double *t, size_t n)
{
	double a;
	double b;
	double c;
	size_t i;
	ssize_t r = 0;

	if (bz == NULL || t == NULL)
		ERR(EFAULT, -1);
	if (o == NULL) {
		for (i = 0; i < n; i++)
			r += t[i] >= 0 && t[i] <= 1;
		return (r);
	}
	a = bz[0] - 2 * bz[2] + bz[4];
	b = 2 * (bz[2] - bz[0]);
	c = bz[0];
	for (i = 0; i < n; i++) {
		if (t[i] < 0 || t[i] > 1)
			continue;
		o[r++] = a * t[i] * t[i] + b * t[i] + c;
	}
	return (r);
}

/* Evaluate for Y at `t' */
ssize_t
bz2y(double *bz, double *o, double *t, size_t n)
{
	double a;
	double b;
	double c;
	size_t i;
	ssize_t r = 0;

	if (bz == NULL || t == NULL)
		ERR(EFAULT, -1);
	if (o == NULL) {
		for (i = 0; i < n; i++)
			r += t[i] >= 0 && t[i] <= 1;
		return (r);
	}
	a = bz[1] - 2 * bz[3] + bz[5];
	b = 2 * (bz[3] - bz[1]);
	c = bz[1];
	for (i = 0; i < n; i++) {
		if (t[i] < 0 || t[i] > 1)
			continue;
		o[r++] = a * t[i] * t[i] + b * t[i] + c;
	}
	return (r);
}

/* Evaluate for X at `t' */
ssize_t
bz3x(double *bz, double *o, double *t, size_t n)
{
	double a;
	double b;
	double c;
	double d;
	size_t i;
	ssize_t r = 0;

	if (bz == NULL || t == NULL)
		ERR(EFAULT, -1);
	if (o == NULL) {
		for (i = 0; i < n; i++)
			r += t[i] >= 0 && t[i] <= 1;
		return (r);
	}
	a = bz[6] - 3 * bz[4] + 3 * bz[2] - bz[0];
	b = 3 * (bz[4] - 2 * bz[2] + bz[0]);
	c = 3 * (bz[2] - bz[0]);
	d = bz[0];
	for (i = 0; i < n; i++) {
		if (t[i] < 0 || t[i] > 1)
			continue;
		o[r++] = a * t[i] * t[i] * t[i] + b * t[i] * t[i] +
		    c * t[i] + d;
	}
	return (r);
}

/* Evaluate for Y at `t' */
ssize_t
bz3y(double *bz, double *o, double *t, size_t n)
{
	double a;
	double b;
	double c;
	double d;
	size_t i;
	ssize_t r = 0;

	if (bz == NULL || t == NULL)
		ERR(EFAULT, -1);
	if (o == NULL) {
		for (i = 0; i < n; i++)
			r += t[i] >= 0 && t[i] <= 1;
		return (r);
	}
	a = bz[7] - 3 * bz[5] + 3 * bz[3] - bz[1];
	b = 3 * (bz[5] - 2 * bz[3] + bz[1]);
	c = 3 * (bz[3] - bz[1]);
	d = bz[1];
	for (i = 0; i < n; i++) {
		if (t[i] < 0 || t[i] > 1)
			continue;
		o[r++] = a * t[i] * t[i] * t[i] + b * t[i] * t[i] +
		    c * t[i] + d;
	}
	return (r);
}

/* Evaluate the X derivitive at `t' */
ssize_t
bz2xp(double *bz, double *o, double *t, size_t n)
{
	double a;
	double b;
	ssize_t i;
	ssize_t r = 0;

	if (bz == NULL || t == NULL)
		ERR(EFAULT, -1);
	if (o == 0) {
		for (i = 0; i < n; i++)
			r += t[i] >= 0 && t[i] <= 1;
		return (r);
	}
	a = 2 * (bz[0] - 2 * bz[2] + bz[4]);
	b = 2 * (bz[2] - bz[0]);
	for (i = 0; i < n; i++) {
		if (t[i] < 0 || t[i] > 1)
			continue;
		o[r++] = a * t[i] + b;
	}
	return (r);
}

/* Evaluate the Y derivitive at `t' */
ssize_t
bz2yp(double *bz, double *o, double *t, size_t n)
{
	double a;
	double b;
	ssize_t i;
	ssize_t r = 0;

	if (bz == NULL || t == NULL)
		ERR(EFAULT, -1);
	if (o == 0) {
		for (i = 0; i < n; i++)
			r += t[i] >= 0 && t[i] <= 1;
		return (r);
	}
	a = 2 * (bz[1] - 2 * bz[3] + bz[5]);
	b = 2 * (bz[3] - bz[1]);
	for (i = 0; i < n; i++) {
		if (t[i] < 0 || t[i] > 1)
			continue;
		o[r++] = a * t[i] + b;
	}
	return (r);
}

/* Evaluate the X derivitive at `t' */
ssize_t
bz3xp(double *bz, double *o, double *t, size_t n)
{
	double a;
	double b;
	double c;
	ssize_t i;
	ssize_t r = 0;

	if (bz == NULL || t == NULL)
		ERR(EFAULT, -1);
	if (o == 0) {
		for (i = 0; i < n; i++)
			r += t[i] >= 0 && t[i] <= 1;
		return (r);
	}
	a = 3 * (bz[6] - 3 * bz[4] + 3 * bz[2] - bz[0]);
	b = 6 * (bz[4] - 2 * bz[2] + bz[0]);
	c = 3 * (bz[2] - bz[0]);
	for (i = 0; i < n; i++) {
		if (t[i] < 0 || t[i] > 1)
			continue;
		o[r++] = a * t[i] * t[i] + b * t[i] + c;
	}
	return (r);
}

/* Evaluate the Y derivitive at `t' */
ssize_t
bz3yp(double *bz, double *o, double *t, size_t n)
{
	double a;
	double b;
	double c;
	ssize_t i;
	ssize_t r = 0;

	if (bz == NULL || t == NULL)
		ERR(EFAULT, -1);
	if (o == 0) {
		for (i = 0; i < n; i++)
			r += t[i] >= 0 && t[i] <= 1;
		return (r);
	}
	a = 3 * (bz[7] - 3 * bz[5] + 3 * bz[3] - bz[1]);
	b = 6 * (bz[5] - 2 * bz[3] + bz[1]);
	c = 3 * (bz[3] - bz[1]);
	for (i = 0; i < n; i++) {
		if (t[i] < 0 || t[i] > 1)
			continue;
		o[r++] = a * t[i] * t[i] + b * t[i] + c;
	}
	return (r);
}

/* Calculate the local maximum/minimum value for X */
ssize_t
bz2xl(double *bz, double *to, double *x, double *y)
{
	double a;
	double b;
	double t;

	if (bz == NULL)
		ERR(EFAULT, -1);
	a = bz[0] - 2 * bz[2] + bz[4];
	b = 2 * (bz[2] - bz[0]);
	if (a == 0)
		return (0);
	t = -b / a / 2;
	if (t < 0 || t > 1)
		return (0);
	if (to != NULL)
		*to = t;
	if (x != NULL)
		*x = a * t * t + b * t + bz[0];
	if (y != NULL)
		(void)bz2y(bz, x, &t, 1);
	return (1);
}

/* Calculate the local maximum/minimum value for Y */
ssize_t
bz2yl(double *bz, double *to, double *x, double *y)
{
	double a;
	double b;
	double t;

	if (bz == NULL)
		ERR(EFAULT, -1);
	a = bz[1] - 2 * bz[3] + bz[1];
	b = 2 * (bz[3] - bz[1]);
	if (a == 0)
		return (0);
	t = -b / a / 2;
	if (t < 0 || t > 1)
		return (0);
	if (to != NULL)
		*to = t;
	if (y != NULL)
		*y = a * t * t + b * t + bz[1];
	if (x != NULL)
		(void)bz2x(bz, x, &t, 1);
	return (1);
}

/* Calculate the X intercept at `y' */
ssize_t
bz1xi(double *bz, double *to, double *x, double y)
{
	double t;
	double a;
	double b;

	if (bz == NULL)
		ERR(EFAULT, -1);
	if ((y < bz[1] && y < bz[3]) ||
	    (y > bz[1] && y > bz[3]))
		return (0);
	a = bz[3] - bz[1];
	b = bz[1];
	if (a == 0)
		return (0);
	t = -b / a;
	if (t < 0 || t > 1)
		return (0);
	if (to != NULL)
		to[0] = t;
	if (x != NULL)
		(void)bz1x(bz, x, &t, 1);
	return (1);
}

/* Calculate the Y intercept at `x' */
ssize_t
bz1yi(double *bz, double *to, double x, double *y)
{
	double t;
	double a;
	double b;

	if (bz == NULL)
		ERR(EFAULT, -1);
	if ((x < bz[0] && x < bz[2]) ||
	    (x > bz[0] && x > bz[2]))
		return (0);
	a = bz[2] - bz[0];
	b = bz[0];
	if (a == 0)
		return (0);
	t = -b / a;
	if (t < 0 || t > 1)
		return (0);
	if (to != NULL)
		*to = t;
	if (y != NULL)
		(void)bz1y(bz, y, &t, 1);
	return (1);
}

/* Calculate the X intercept at `y' */
ssize_t
bz2xi(double *bz, double *to, double *x, double y, int tflag)
{
	double t[2];
	ssize_t r;

	if (bz == NULL)
		ERR(EFAULT, -1);
	if ((y < bz[1] && y < bz[3] && y < bz[5]) ||
	    (y > bz[1] && y > bz[3] && y > bz[5]))
		return (0);
	r = quadratic(bz[1] - 2 * bz[3] + bz[5], 2 * (bz[3] - bz[1]),
	    bz[1] - y, t, tflag);
	switch (r) {
	case 1:
		if (t[0] < 0 || t[0] > 1)
			r--;
		break;
	case 2:
		if (t[1] < 0 || t[1] > 1)
			r--;
		if (t[0] < 0 || t[0] > 1) {
			t[0] = t[1];
			r--;
		}
		break;
	}
	if (r > 0) {
		if (to != NULL)
			(void)memcpy(to, t, r * sizeof(double));
		if (x != NULL)
			(void)bz2x(bz, x, t, 2);
	}
	return (r);
}

/* Calculate the Y intercept at `x' */
ssize_t
bz2yi(double *bz, double *to, double x, double *y, int tflag)
{
	double t[2];
	ssize_t r;

	if (bz == NULL)
		ERR(EFAULT, -1);
	if ((x < bz[0] && x < bz[2] && x < bz[4]) ||
	    (x > bz[0] && x > bz[2] && x > bz[4]))
		return (0);
	r = quadratic(bz[0] - 2 * bz[2] + bz[4], 2 * (bz[2] - bz[0]),
	    bz[0] - x, t, tflag);
	switch (r) {
	case 1:
		if (t[0] < 0 || t[0] > 1)
			r--;
		break;
	case 2:
		if (t[1] < 0 || t[1] > 1)
			r--;
		if (t[0] < 0 || t[0] > 1) {
			t[0] = t[1];
			r--;
		}
		break;
	}
	if (r > 0) {
		if (to != NULL)
			(void)memcpy(to, t, r * sizeof(double));
		if (y != NULL)
			(void)bz2y(bz, y, t, 2);
	}
	return (r);
}

/* Cut the curve at the specifiend `t' values */
ssize_t
bz2c(double *bz, double ta, double tb, double *o)
{
	double tta;
	double ttb;
	double ttd;
	double ax;
	double ay;
	double bx;
	double by;
	double td;

	if (bz == NULL)
		ERR(EFAULT, -1);
	if (o == NULL)
		return (6);
	if (ta == 0 && tb == 1) {
		(void)memcpy(o, bz, 6 * sizeof(double));
		return (6);
	}
	ax = bz[0] - 2 * bz[2] + bz[4];
	ay = bz[1] - 2 * bz[3] + bz[5];
	bx = 2 * (bz[2] - bz[0]);
	by = 2 * (bz[3] - bz[1]);
	tta = ta * ta;
	ttb = tb * tb;
	if (ta == tb) {
		o[0] = ax * tta + bx * ta + bz[0];
		o[1] = ay * tta + by * ta + bz[1];
		(void)memcpy(o + 2, o, 2 * sizeof(double));
		(void)memcpy(o + 4, o, 2 * sizeof(double));
		return (6);
	}
	td = (tb + ta) / 2;
	ttd = td * td;
	o[0] = ax * tta + bx * ta + bz[0];
	o[1] = ay * tta + by * ta + bz[1];
	o[2] = ax * ttd + bx * td + bz[0];
	o[3] = ay * ttd + by * td + bz[1];
	o[4] = ax * ttb + bx * tb + bz[0];
	o[5] = ay * ttb + by * tb + bz[1];
	o[2] = 2 * o[2] - o[0] / 2 - o[4] / 2;
	o[3] = 2 * o[3] - o[1] / 2 - o[5] / 2;
	return (6);
}

/*
 * Calculate the area of the curve segment.
 *
 * The equation for the area was obtained by the sum of the area of an
 * an infinite series of infinitesimal triangles covering the area of
 * the curve.
 */
ssize_t
bz2a(double *bz, double *a)
{

	if (bz == NULL)
		ERR(EFAULT, -1);
	if (a != NULL)
		*a = (-bz[0] * bz[3] + bz[0] * bz[5] + bz[2] * bz[1] -
		    bz[2] * bz[5] - bz[4] * bz[1] + bz[4] * bz[3]) / 3;
	return (1);
}

/*
 * See above for the method used.
 * Note: this has not been tested
 */
ssize_t
bz3a(double *bz, double *a)
{

	if (bz == NULL)
		ERR(EFAULT, -1);
	if (a != NULL)
		*a = 3 * (2 * bz[0] * bz[3] + bz[0] * bz[5] -
		    3 * bz[0] * bz[7] - 2 * bz[2] * bz[1] +
		    bz[2] * bz[5] + bz[2] * bz[7] - bz[4] * bz[1] -
		    bz[4] * bz[3] + 2 * bz[4] * bz[7] +
		    3 * bz[6] * bz[1] - bz[6] * bz[3] -
		    2 * bz[6] * bz[5]) / 20;
	return (1);
}

/* Calculate the area of quadratic curves */
ssize_t
bz2aa(double *bz, size_t count, double *ao)
{
	double a;
	double aa;
	double t;
	size_t i;
	size_t j;

	if (bz == NULL)
		ERR(EFAULT, -1);
	if (ao == NULL)
		return (1);
	if (count < 6) {
		*ao = 0;
		return (1);
	}
	for (i = 0, j = 0, a = 0, aa = 0; i + 6 <= count; i += 6) {
		(void)bz2a(bz + i, &t);
		a += t;
		aa += (bz[i + 0] * bz[i + 5] -
		    bz[i + 4] * bz[i + 1]) / 2;
	}
	*ao = aa + a;
	return (1);
}

ssize_t
bz2split(double *bz, size_t bc, double **od, size_t *oc, size_t oi,
    double v, int fl)
{
	double d[6];
	double h[6];
	double td[2];
	double mod;
	double t;
	size_t imod;
	size_t tc;
	size_t i;
	int side;
	int f;

	if (bz == NULL)
		ERR(EFAULT, -1);
	if (bc < 6)
		return (0);
	mod = (fl & BZ_LOW) ? -1 : 1;
	imod = (fl & BZ_Y) ? 1 : 0;
	/* Find the initial side of `v'*/
	side = (bz[imod] - v) > 0;
	/* Follow the curve looking for crossings */
	for (i = 0, f = 0; i + 6 <= bc; i += 6) {
		d[0] = (bz[i + 0 + imod] - v) * mod;
		d[1] = (bz[i + 2 + imod] - v) * mod;
		d[2] = (bz[i + 4 + imod] - v) * mod;
		if ((d[0] < 0 && d[1] < 0 && d[2] < 0) ||
		    (d[0] == 0 && d[1] < 0 && d[2] == 0)) {
			continue;
		}
		if (d[0] > 0 && d[1] > 0 && d[2] > 0) {
			/* This segment is OK */
			if (od != NULL && addvalues(od, oc, oi, bz + i,
			    6) < 0)
				return (-1);
			oi += 6;
			continue;
		}
		/* This segment crosses the boundry */
		if (imod)
			tc = bz2xi(bz + i, td, NULL, v, 0);
		else
			tc = bz2yi(bz + i, td, v, NULL, 0);
		switch (tc) {
		case 0:
			if (d[0] < 0 || d[2] < 0)
				break;
			if (od != NULL && addvalues(od, oc, oi,
			    bz + i, 6) < 0)
				return (-1);
			oi += 6;
			break;
		case 1:
			if (d[0] > d[2]) {
				/* Exiting the area */
				if (bz2c(bz + i, 0, td[0], d) != 6)
					break;
				if (d[0] != d[2] || d[2] != d[4] ||
				    d[1] != d[3] || d[3] != d[5]) {
					if (od != NULL && addvalues(od,
					    oc, oi, d, 6) < 0)
						return (-1);
					oi += 6;
				}
				if (f && td[0] != 0) {
					h[4] = d[4];
					h[5] = d[5];
					h[2] = (h[0] + h[4]) / 2;
					h[3] = (h[1] + h[5]) / 2;
					t = h[0];
					h[0] = h[4];
					h[4] = t;
					t = h[1];
					h[1] = h[5];
					h[5] = t;
					if (h[0] != h[4] ||
					    h[1] != h[5]) {
						if (od != NULL &&
						    addvalues(od, oc,
						    oi, h, 6) < 0)
							return (-1);
						oi += 6;
					}
				} else if (td[0] != 0) {
					h[0] = d[4];
					h[1] = d[5];
				}
				if (td[0] != 0)
					f = !f;
			} else {
				/* Entering the area */
				if (bz2c(bz + i, td[0], 1, d) != 6)
					break;
				if (f && td[0] != 1) {
					h[4] = d[0];
					h[5] = d[1];
					h[2] = (h[0] + h[4]) / 2;
					h[3] = (h[1] + h[5]) / 2;
					if (h[0] != h[4] ||
					    h[1] != h[5]) {
						if (od != NULL &&
						    addvalues(od, oc,
						    oi, h, 6) < 0)
							return (-1);
						oi += 6;
					}
				} else if (td[0] != 1) {
					h[0] = d[0];
					h[1] = d[1];
				}
				if (td[0] != 1)
					f = !f;
				if (d[0] != d[2] || d[2] != d[4] ||
				    d[1] != d[3] || d[3] != d[5]) {
					if (od != NULL && addvalues(od,
					    oc, oi, d, 6) < 0)
						return (-1);
					oi += 6;
				}
			}
			break;
		case 2:
			if (td[0] > td[1]) {
				d[3] = td[0];
				td[0] = td[1];
				td[1] = d[3];
			}
			if (d[1] > 0) {
				/* Pops in & out */
				if (bz2c(bz + i, td[0], td[1],
				    d) != 6)
					break;
				if (f) {
					h[4] = d[0];
					h[5] = d[1];
					h[2] = (h[0] + h[4]) / 2;
					h[3] = (h[1] + h[5]) / 2;
					if (h[0] != h[4] ||
					    h[1] != h[5]) {
						if (od != NULL &&
						    addvalues(od, oc,
						    oi, h, 6) < 0)
							return (-1);
						oi += 6;
					}
				} else {
					h[0] = d[0];
					h[1] = d[1];
				}
				f = !f;
				if (d[0] != d[2] || d[2] != d[4] ||
				    d[1] != d[3] || d[3] != d[5]) {
					if (od != NULL && addvalues(od,
					    oc, oi, d, 6) < 0)
						return (-1);
					oi += 6;
				}
				if (f) {
					h[4] = d[4];
					h[5] = d[5];
					h[2] = (h[0] + h[4]) / 2;
					h[3] = (h[1] + h[5]) / 2;
					t = h[0];
					h[0] = h[4];
					h[4] = t;
					t = h[1];
					h[1] = h[5];
					h[5] = t;
					if (h[0] != h[4] ||
					    h[1] != h[5]) {
						if (od != NULL &&
						    addvalues(od, oc,
						    oi, h, 6) < 0)
							return (-1);
						oi += 6;
					}
				} else {
					h[0] = d[4];
					h[1] = d[5];
				}
				f = !f;
			} else {
				/* Pops out & in */
				if (bz2c(bz + i, 0, td[0], d) != 6)
					break;
				if (d[0] != d[2] || d[2] != d[4] ||
				    d[1] != d[3] || d[3] != d[5]) {
					if (od != NULL && addvalues(od,
					    oc, oi, d, 6) < 0)
						return (-1);
					oi += 6;
				}
				if (f && td[0] != 0) {
					h[4] = d[4];
					h[5] = d[5];
					h[2] = (h[0] + h[4]) / 2;
					h[3] = (h[1] + h[5]) / 2;
					t = h[0];
					h[0] = h[4];
					h[4] = t;
					t = h[1];
					h[1] = h[5];
					h[5] = t;
					if (h[0] != h[4] ||
					    h[1] != h[5]) {
						if (od != NULL &&
						    addvalues(od, oc,
						    oi, h, 6) < 0)
							return (-1);
						oi += 6;
					}
				} else if (td[0] != 0) {
					h[0] = d[4];
					h[1] = d[5];
				}
				if (td[0] != 0)
					f = !f;
				if (bz2c(bz + i, td[1], 1, d) != 6)
					break;
				if (f && td[1] != 1) {
					h[4] = d[0];
					h[5] = d[1];
					h[2] = (h[0] + h[4]) / 2;
					h[3] = (h[1] + h[5]) / 2;
					if (h[0] != h[4] ||
					    h[1] != h[5]) {
						if (od != NULL &&
						    addvalues(od, oc,
						    oi, h, 6) < 0)
							return (-1);
						oi += 6;
					}
				} else if (td[1] != 1){
					h[0] = d[0];
					h[1] = d[1];
				}
				if (td[1] != 1)
					f = !f;
				if (d[0] != d[2] || d[2] != d[4] ||
				    d[1] != d[3] || d[3] != d[5]) {
					if (od != NULL && addvalues(od,
					    oc, oi, d, 6) < 0)
						return (-1);
					oi += 6;
				}
			}
			break;
		}
	}
	return (oi);
}

int
addvalues(double **o, size_t *oc, size_t oi, double *d, size_t dc)
{
	size_t c;
	void *p;

	if (o == NULL)
		ERR(EFAULT, -1);
	if (d == NULL || dc == 0)
		return (0);
	if (oc != NULL) {
		c = *oc;
		if (c < oi)
			ERR(EFAULT, -1);
	} else
		c = oi;
	if (c - oi < dc) {
		if ((p = reallocarray(*o, c + dc,
		    sizeof(double))) == NULL)
			return (-1);
		*o = p;
		if (oc != NULL)
			*oc = c + dc;
	}
	(void)memcpy(*o + oi, d, dc * sizeof(double));
	return (0);
}

ssize_t
quadratic(double a, double b, double c, double *o, int tflag)
{
	double s;

	if (o == NULL)
		ERR(EFAULT, -1);
	if (a == 0) {
		if (b == 0)
			return (0);
		o[0] = -c / b;
		return (1);
	}
	s = b * b - 4 * a * c;
	if (s < 0)
		return (0);
	if (s == 0) {
		if (!tflag)
			return (0);
		o[0] = b / a;
		return (1);
	}
	s = -b + (b < 0 ? 1 : -1) * sqrt(s);
	o[0] = s / a / 2;
	o[1] = c / s * 2;
	return (2);
}
